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We present a systematic analysis of two-pion interferometry in Au+Au collisions at ^/snn = 200 
GeV using the STAR detector at RHIC. We extract the HBT radii and study their multiplicity, 
transverse momentum, and azimuthal angle dependence. The Gaussianess of the correlation func- 
tion is studied. Estimates of the geometrical and dynamical structure of the freeze-out source are 
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extracted by fits with blast wave parameterizations. The expansion of the source and its relation 
with the initial energy density distribution is studied. 



1. INTRODUCTION 

In the late 1950's two particle intensity inter- 
ferometry (HBT) was proposed and developed 
by the astronomers Hanbury-Brown and Twiss 
to measure the angular size of distant stars Q . 
In 1960, Goldhaber et al. applied this tech- 
nique to particle physics to study the angular 
distribution of identical pion pairs in pp anni- 
hilations 0- They observed an enhancement 
of pairs at small relative momenta that was ex- 
plained in terms of the symmetrization of the 
two-pion wave function. 

In ultra-relativistic heavy ion collisions, 
where a quark gluon plasma (QGP) is expected 
to be formed, HBT is a useful tool to study the 
space-time geometry of the particle-emitting 
source [1 El It also contains dynamical infor- 
mation that can be explored by studying the 
transverse momentum dependence of the ap- 
parent source size 0.0]- In non-central colli- 
sions, information on the anisotropic shape of 
the pion-emitting region at kinetic freeze-out 
can be extracted by measuring two-pion corre- 
lation functions as a function of the emission 
angle with respect the reaction plane, see, for 
example, 0,H,|^|. 

Experimentally, two-particle correlations are 
studied by constructing the correlation function 
as H: 



C(q) 



B(q)' 



(1) 



Here A(q) is the pair distribution in momentum 
difference q = p\ — p\ for pairs of particles from 
the same event, and B(q) is the corresponding 
distribution for pairs of particles from differ- 
ent events. To good approximation this ratio 
is sensitive to the spatial extent of the emit- 
ting region and insensitive to the single particle 
momentum distribution, acceptance, and effi- 
ciency effects 0. 

At the Relativistic Heavy Ion Collider 
(RHIC), identical-pion HBT studies at ^snn 
= 130 GeV 0, lead to an apparent 
source size qualitatively similar to measure- 
ments at lower energies. In contrast to pre- 
dictions of larger sources based on QGP for- 
mation [l2l Il3j |. no long emission duration is 
seen. The extracted parameters do not agree 
with predictions of hydrodynamic models that, 
on the other hand, describe reasonably well 
the momentum-space structure of the emitting 
source and elliptic flow [3 . This "HBT puzzle" 



could be related to the fact that the extracted 
timescales are smaller than those predicted by 
the hydrodynamical model [l3j|. More sophis- 
ticated approaches such as 3D hydrodynami- 
cal calculations 0] or multi-stage models 0] 
also cannot describe simultaneously the geom- 
etry and the dynamic of the system 0] . 



Further detail may be obtained from non- 
central collisions, where the initial anisotropic 
collision geometry has an almond shape with its 
longer axis perpendicular to the reaction plane. 
This generates greater transverse pressure gra- 
dients in the reaction plane than perpendicular 
to it. This leads to preferential in-plane expan- 
sion [lil Il9l I20I l2lf which diminishes the ini- 
tial anisotropy as the source evolves. Thus, the 
source shape at freeze-out should be sensitive 



to the evolution of the pressure and th e sy stem 
lifetime. Hydrodynamic calculations [22| pre- 
dict that the source may still be out-of-plane 
extended after hydrodynamic evolution. How- 
ever a subsequent rescattering phase tends to 
make the source in-plane extended [2^ . There- 
fore, the experimental freeze-out source shape 
might discriminate between different scenarios 
of the system's evolution. 



In this paper we present results of our sys- 
tematic studies of two-pion HBT correlations 
in Au+Au collisions at ^/snn — 200 GeV 
measured in the STAR (Solenoidal Tracker at 
RHIC) detector at RHIC. We describe the anal- 
ysis procedure in detail and discuss several is- 
sues with importance to HBT such as different 
ways of taking the final state Coulomb inter- 
action into account, and the Gaussianess of the 
measured correlation function. The paper is or- 
ganized as follows: Section introduces the ex- 
perimental setup as well as the event, particle 
and pair selections. In section [3] the analy- 
sis method is presented. In section system- 
atic results are shown. We discuss these re- 
sults in section [13 where the centrality depen- 
dence of the transverse mass rnx dependence 
of the HBT parameters is investigated, the ex- 
tracted parameters from a fit to a blast wave 
parametrization are discussed in detail and the 
expansion of the source is studied. We summa- 
rize and conclude in sectional Extended details 
about the analysis method, the results and the 
discussion can be found in reference 241. 
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2. EXPERIMENTAL SETUP, EVENT 
AND PARTICLE SELECTION 

2.1. STAR detector 

The STAR detector is an azimuthally sym- 
metric, large acceptance, solenoidal detector. 
The subsystems relevant for this analysis are 
a large Time Projection Chamber (TPC) lo- 
cated in a 0.5 Tesla solenoidal magnet, two zero- 
degree calorimeters (ZDCs) that detect spec- 
tator neutrons from the collision, and a cen- 
tral trigger barrel (CTB) that measures charged 
particle multiplicity. The latter two subsystems 
were used for online triggering only. 

The TPC [13 is the primary STAR detec- 
tor and the only detector used for the event 
reconstruction of the analysis presented here. 
It is 4.2 m long and covers the pseudorapidity 
region \rj\ < 1.8 with full azimuthal coverage 
(— 7r < cf> < 7r). It is a gas chamber, with inner 
and outer radii of 50 and 200 cm respectively, 
in a uniform electric field. The particles pass- 
ing through the gas release secondary electrons 
that drift to the readout end caps at both ends 
of the chamber. The readout system is based on 
multiwirc proportional chambers, with readout 
pads. There are 45 pad- rows between the inner 
and the outer radii of the TPC. The induced 
charge from the electrons is shared over several 
adjacent pads, so the original track position is 
reconstructed to ~ 500 /im precision. 

The STAR trigger detectors are the CTB and 
the ZDCs. In this analysis two trigger settings 
were used. Hadronic minimum-bias that re- 
quires a signal above threshold in both ZDCs, 
and hadronic central that requires low ZDC sig- 
nal and high CTB signal. 



2.2. Event selection and binning 

For this analysis, we selected events with a 
collision vertex position within ± 25 cm mea- 
sured along the beam axis from the center of 
the TPC. This event selection was applied to 
all data sets discussed here. 

We further binned events by centrality, where 
the centrality was characterized according to 
the measured multiplicity of charged hadrons 
with pseudorapidity (|ry| < 0.5), and here we 
present results as a function of centrality bins. 
The six centrality bins correspond to 0-5%, 5- 
10%, 10-20%, 20-30%, 30-50%, and 50-80% of 
the total hadronic cross-section. A hadronic- 
central triggered data set of 1 million events 
was used only for the first bin. The other five 



bins are from a minimum-bias triggered data 
set of 1.7 million events. 

Within each centrality bin, in order to form 
the background pairs for the correlation func- 
tion (see section 13.10 . we only mixed "simi- 
lar" events. In this analysis, "similar" events 
have primary vertex relative z position within 
5 cm, multiplicities within the same centrality 
bin described above, and, for the azimuthally- 
sensitive analysis, estimated reaction plane ori- 
entations within 20°. 



2.3. Particle selection 

We selected tracks in the rapidity region \y\ < 
0.5. Particle identification was done by corre- 
lating the specific ionization of the particles in 
the gas of TPC with their measured momen- 
tum. The energy lost by a particle as it trav- 
els through a gas depends on the velocity j3 at 
which it travels and it is described by the Bcthe- 
Bloch formula 26]. 

For a given momentum, each particle mass 
will have a different velocity and a different 
dE/dx as it goes through the gas of the TPC. 
For this analysis pions were selected by requir- 
ing the specific ionization to be within 2 stan- 
dard deviations (experimentally determined as 
a function of the particle momentum and event 
multiplicity) from the Bcthc-Bloch value for pi- 
ons. To help remove kaons that could satisfy 
this condition, particles were also required to 
be farther than 2 standard deviations from the 
value for kaons. There is a small contamina- 
tion of electrons in the low momentum region, 
p < 400 MeV/c. Its effect was studied with 
different cuts and found to be unimportant. 

To reduce contributions from non-primary 
(decay) pions, we applied a cut of 3 cm to each 
track on the distance of closest approach of the 
extrapolated track to the primary vertex. 

In our previous HBT analysis [Toj |. tracks 
were divided in different bins according to their 
transverse momentum, px, and only particles 
within a given bin were used to form correlation 
functions. In this analysis no such px binning 
was applied. The px range was set by limita- 
tions in the reconstruction of pions in the TPC, 
by the fact that we remove the kaon band and 
by the momentum pair cut described below and 
only tracks with 150 < px < 800 MeV/c were 
accepted. 
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a) 



FIG. 1: Distribution of same number of hits in two 
tracks for four possible cases. Closed circles are hits 
assigned to one track, open circles are assigned to 
the other, a) SL = -0.5 (clearly two tracks) b) SL 
= 1 (possible split track) c) SL = 1 (possible split 
track) d) SL = 0.08 (likely two tracks). 



2.4. Pair cuts 

In this section we describe pair cuts and bin- 
ning. The first two cuts discussed are intended 
to remove the effects of two track reconstruc- 
tion defects that are important to HBT: split 
tracks (one single particle reconstructed as two 
tracks) and merged tracks (two particles with 
similar momenta reconstructed as one track). 



2.4-1- Split tracks 

Track splitting causes an enhancement of 
pairs at low relative momentum q. This false 
enhancement is created by single tracks recon- 
structed as two with similar momenta. In order 
to remove split tracks we compare the location 
of the hits for each track in the pair along the 
pad-rows in the TPC and assign a quantity to 
each pair, called Splitting Level (SL), calculated 
as follows: 



SL 



where Si 



(2) 



Nhitsi + Nhits 2 

{+1 one track leaves a hit on pad-row 
— 1 both tracks leave a hit on pad-row 
neither track leaves a hit on pad-row, 

where i is the pad-row number, and Nhitsi and 
NhitS2 are the total number of hits associated 
to each track in the pair. If only one track has 
a hit in a pad-row +1 is added to the running 
quantity, if both tracks have a hit in the same 
pad-row, a sign of separate tracks, — 1 is added 
to this quantity. After the sum is done, it is 
divided by the sum of hits in both tracks, this 
normalizes SL to a value between -0.5 (both 
tracks have hits in exactly same pad-rows) and 
1.0 (tracks do not have any hit in same pad- 
row) . Figure ^ shows four different cases for 
the same number of total hits: in case a) two 
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FIG. 2: ID correlation function for different values 
of SL (anti-splitting cut). The cut applied in the 
analysis is SL < 0.6. The horizontal lines indicate 
the bin width. 



different tracks with SL = -0.5, in b) and c) 
two different cases of possible split tracks with 
SL = 1, and in d) two different tracks with SL 
= 0.08. 

We required every pair to have SL smaller 
than a certain value. This value was determined 
from the 1-dimensional correlation functions as 
a function of the relative momentum of the pair 
<7inv, for different values of SL; some of them are 
shown in Fig. [21 The relative momentum of the 
pair is defined as q- ulv — y (q ) 2 — \q\ 2 where q° 
and q are the components of the four- vector mo- 
mentum difference. We observe that when mak- 
ing this cut more restrictive (reducing the max- 
imum allowed value for SL) the enhancement is 
reduced until we reach SL = 0.6 when the cor- 
relation function becomes stable and does not 
change for lower values of SL. Therefore, all the 
pairs entering the correlation functions were re- 
quired to have SL < 0.6. Cutting at this value 
is also supported by simulation studies. While 
naturally track splitting can only give rise to 
false pairs, our SL cut also removes some real 
pairs which happen to satisfy the cut. There- 
fore, we apply the SL cut to both, "real" and 
"mixed" pairs, numerator and denominator of 
C(q), Eq. Q. 



2-4-2. Merged tracks 

Once we have removed split tracks we can 
study the effects of two particles reconstructed 
as one track. These merged tracks cause a re- 
duction of pairs at low relative momentum since 
the particles that have higher probability of be- 
ing merged are those with similar momenta. To 
eliminate the effect of track merging, we re- 
quired that all pairs entering numerator and 
denominator of the correlation function had a 
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FIG. 3: ID correlation function for different values 
of the maximum fraction of merged hits allowed. 
Cut applied in the analysis is fraction of merged 
hits (FMH) < 10%. The horizontal lines indicate 
the bin width. 



fraction of merged hits no larger than 10%. 
Two hits are considered merged if the probabil- 
ity of separating them is less than 99%. From 
simulation and data studies, this minimum sep- 
aration was determined to be 5 mm. By apply- 
ing this cut to "real" and "mixed" pairs, we 
introduce in the denominator the effect that 
merged tracks have in the numerator: a reduc- 
tion of low q pairs. Note that tracks from differ- 
ent events will originate from primary vertices 
at different positions along the beam direction. 
Thus, even two tracks with identical momenta, 
which would surely be merged if they originated 
from the same event, may not be considered as 
a merged track when they formed a "mixed" 
pair if we would not account for the different 
primary vertex position. Our procedure is to 
calculate (using a helix model) the pad-row hit 
positions of each track, assuming that the track 
originated at the center of the TPC (24|. These 
calculated hit positions are used in the merg- 
ing cut procedure described above. By apply- 
ing this cut to the numerator, we would remove 
"real" pairs that satisfy the cut. This would 
reduce the HBT fit parameters for a correla- 
tion function that is not completely Gaussian, 
and needs to be taken into account as will be 
described in the next section. 

To determine the maximum fraction of 
merged hits allowed we proceed as we did for 
the anti-splitting cut. Figure shows the 1- 
dimcnsional correlation functions as a function 
of gi n v, for different values of the maximum frac- 
tion of merged hits allowed. By requiring the 
fraction of merged hits to be less than 10% for 
every pair entering the correlation function, the 
effect of merged tracks in the correlation func- 
tion was almost completely removed as will be 
discussed in section |3] 



2.4-3. kr cut and pair binning 



As already mentioned, no explicit px cut was 
applied to single tracks beyond the require- 
ment for clean PID. However, in addition to 
the two cuts already described, pairs were re- 
quired to have an average transverse momen- 
tum (kr = (\pit\ + |p2t|)/2) between 150 and 
600 MeV/c. No difference was observed be- 
tween the extracted HBT parameters when ap- 
plying equivalent px or kx cuts. However statis- 
tics improved when using the latter cut, as two 
pions from different pt bins will be used in a 
fcr-cut analysis, but not in a px-cut analysis. 

Pairs were then binned by kx in 4 bins 
that correspond to [150,250] MeV/c, [250,350] 
McV/c, [350,450] MeV/c, and [450,600] MeV/c. 
Here the results are presented as a function of 
the average kx (or mj = y/k T + m*) in each 
of those bins. 

In addition, in the azimuthally-sensitive HBT 
analysis, to observe the particle source from 
a series of angles, pairs were also binned ac- 
cording to the angle <& = </> pa ir — where 
0pair is the azimuthal angle of the pair trans- 
verse momentum kx and 'J 2 is the second-order 
event plane azimuthal angle. The first-order 
event plane angle is not reconstructible with the 
STAR detector configuration for this analysis. 
Because we use the second-order reaction plane, 
$ is only defined in the range [0,7r]. 



ANALYSIS METHOD 



3.1. Construction of correlation function 



The two-particle correlation function be- 
tween identical bosons with momenta p\ and 
p-2 is defined in Eq. (Q. As already mentioned, 
A(q) is the measured distribution of the mo- 
mentum difference for pairs of particles from 
the same event and B(q) is obtained by mix- 
ing particles in separate events |27| and repre- 
sents the product of single particle probabili- 
ties. Each particle in one event is mixed with 
all the particles in a collection of events which 
in our case consists of 20 events. As discussed 
before, events in a given collection have pri- 
mary vertex z position within 5 cm, multiplici- 
ties within 5 to 30% of each other, and, for the 
azimuthally-sensitive analysis, estimated reac- 
tion plane orientations within 20°. 
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3.2. Pratt-Bertsch parametrization 

In order to probe length scales differentially 
in beam and transverse directions, the rela- 
tive momentum q is usually decomposed in the 
Pratt-Bertsch (or "out-side-long") convention 
j23,|2^,|23|. In this parametrization the relative 
momentum vector of the pair q is decomposed 
into a longitudinal direction along the beam 
axis, qi, an outward direction parallel to the 
pair transverse momentum, q 0l and a sideward 
direction perpendicular to those two, q s . 

We choose as the reference frame, the lon- 
gitudinal comoving system (LCMS) frame of 
the pair, in which the longitudinal component 
of the pair velocity vanishes. At midrapidity, 
in the LCMS frame, and with knowledge of 
the second-order but not the first-order reac- 
tion plane, the correlation function is usually 
parameterized by a 3-dimensional Gaussian in 
the relative momentum components as [9(: 

C(q) = 1 + ^-tf*S-92J$-tf*?-2«.«.*L ( 3 ) 



Rl(k T ,$) = Rl (k T ) + 

2 J2 <n(McOs(7l$) = 0,8,1) 

n=2,4,6--- 

2 <„(Msin(n$) (/z = os), (4) 

n=2,4,6--- 

where R^^Q^t) are the n th order Fourier coeffi- 
cients for the n radius. These coefficients, that 
are <& independent, can be calculated as: 

2 f (i?2(fc T; $)cos(n$)> (p = o,s,l) 

u ^t) - j (^2 ( fcT; $ j sin ( n $)) = os ). 

(5) 

As we will show, the th order Fourier co- 
efficients correspond to the extracted HBT 
radii in an azimuthally integrated analysis. In 
this analysis we found that Fourier coefficients 
above 2 nd order are consistent with 0. 



For an azimuthally integrated analysis, the cor- 
relation function is symmetric under q s — > —q s 
and R% s = 0. 

In principle, the possibility that the emission 
of particles is neither perfectly chaotic nor com- 
pletely coherent can be taken into account by 
adding the parameter A to the correlation func- 
tion, which, in general, depends on kx- This 
A parameter should be unity for a fully chaotic 
source and smaller than unity for a source with 
partially coherent particle emission. In the 
analysis presented here we have assumed com- 
pletely chaotic emission [3l| and attribute the 
deviations from C(q = 0, k) = 2 to contribution 
from pions coming from long-lived resonances, 
and misidentified particles, such as electrons. 

While for the azimuthally integrated analysis 
the sign of the q components is arbitrary, in the 
azimuthally-sensitive analysis, the sign of R^ s 
is important because it tells us the azimuthal 
direction of the emitted particles, so the signs 
of q and q s are kept and particles in every pair 
are ordered such that qi > 0. 

References [13, H3| give a detailed description 
of the relation between the HBT radius param- 
eters (R% R 2 a , Rf and R^ s ) and the space-time 
geometry of the final freeze-out stage. 



3.3. Fourier components 

For a boost-invariant system, the $ depen- 
dence of the HBT radii of Eq. © are @: 



3.4. Coulomb interaction and fitting 
procedures 

Equation JSJ applies only if the sole cause of 
correlation is quantum statistics and the corre- 
lation function is Gaussian. We come to this 
second point in section |4~T1 In our case, signifi- 
cant Coulomb effects must also be accounted for 
(stronginteractions are within reasonable limit 
here [33]). This Coulomb interaction between 
pairs, repulsive for like-sign particles, causes a 
reduction in the number of real pairs at low q 
reducing the experimental correlation function 
as seen in Fig. ^ 

3.4-1- Standard procedure 

Three different procedures can be applied 
in order to take this interaction into account. 
O ne pro cedure that was used in our analysis at 
\JsNN — 130 GeV 01 as well as by previous 
experiments, consists of fitting the correlation 
function to 

C{q ,q s ,qi) = = K coul (q inv )x (6) 

(1 + \ e -£ R o-<£ R2 '-<ii R i- 2 <i°<i- R lA, 

normalized to unity at large q, where -Kcoui is 
the squared Coulomb wave-function integrated 
over the whole source, which in our case is a 
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FIG. 4: Projections of the 3 dimensional correla- 
tion function and corresponding fits for negative 
pions from the 0-5% most central events and kr 
= [150,250] MeV/c according to the standard and 
Bowler-Sinyukov procedures. 



spherical Gaussian source of 5 fm radius. The 
effect on the final results of changing the radius 
of the spherical Gaussian source to calculate 
-Kcoui was studied and found to be within rea- 
sonable limits. Traditionally, Eq. © has been 
expressed as 



C'(q ,q s ,qi) 



= 1 + Ae" 



A{q) 



B(q)K c . 

■<llK-1 2 3 R2 S- 



iul(<7ii 

i 2 R 2 



-2q„q s R 2 OB 



(7) 



and this new correlation function was called the 
Coulomb corrected correlation function since we 
introduce in the denominator a Coulomb factor 
with which we try to compensate the Coulomb 
interaction in the numerator. We call this stan- 
dard procedure. However, this procedure over- 
corrects the correlation function since it as- 
sumes that all pairs in the background are pri- 
mary pairs and need to be corrected [34[ . 



3.4-2. Dilution procedure 

In a second procedure, inspired by the previ- 
ous procedure and implemented before by the 
E802 collaboration 35], the Coulomb term is 
"diluted" according to the fraction of pairs that 
Coulomb interact as: 

<oul(<Zmv) = 1 + f(K coul (q inv ) - 1), (8) 

where / has a value between (no Coulomb 
weighting) and 1 (standard weight). The corre- 



lation function in this procedure is fitted to: 
A{q) 



C(q ,q s ,qi) 



B{q] 



#coul(<7inv)x (9) 



(1 + \ e -£Rl-fiRl-tiRi- , *<io<isRls\ j 



normalized to unity at large q. We call this the 
dilution procedure. A reasonable assumption is 
to take f — X assuming that A is the fraction 
of primary pions. This increases R by 10-15% 
and has a very small effect on R s and Ri as seen 
in Fig. EI A decreases by 10-15%. 



3.4.3. Bowler-Sinyukov procedure 

An advantage of the previous two techniques 
is that after "correcting" for Coulomb effects, 
one winds up with a correlation function which 
may be fit with a simple Gaussian form. How- 
ever if there exists more than one source of in- 
teraction, it is not valid to "correct" one way. 
For example, it is, in fact, the same pion pairs 
which Coulomb interact as which show quan- 
tum enhancement. This leads to a change in 
the expected form of the correlation function if 
not all particles participate in the interaction 
(i.e., A f 0) HI- If A = 1, all three methods 
are equivalent. 

In this analysis, we have implemented a new 
procedure, first suggested by Bowler [37} and 
Sinyukov et al. [38j , and recently advocated by 
the CERES collaboration (3^, in which only 
pairs with Bose-Einstein interaction are con- 
sidered to Coulomb interact. The correlation 
function in this procedure is fitted to: 



C(q ,q s ,qi) = ^ = (l-X)+ 

A-fr c oul(<7inv)(l + e 



(10) 



normalized to unity at large q, where 
i^coui(9inv) is the same as in the standard pro- 
cedure. The first term on the right-hand side 
of Eq. H10|) accounts for the pairs that do not 
interact and the second term for the pairs that 
(Coulomb and Bose-Einstein) interact. We call 
this Bowler-Sinyukov procedure. It has a simi- 
lar effect on the HBT parameters as the dilution 
procedure as seen in Fig. [5] A similar procedure 
has been recently implemented by the Phobos 
collaboration |64[. In this procedure only pairs 
which are close in the pair center of mass frame 
are considered to Coulomb interact. 

It is worth mentioning that the parameters A 
and and consequently the ratio R /R s , ex- 
tracted using the standard procedure here are 
smaller than the parameters obtained in our 
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FIG. 5: HBT parameters for the three possible 
fitting procedures to the correlation functions de- 
scribed in this paper depending on how Coulomb 
interaction is taken into account from the 0-5% 
most central events. Error bars contain only statis- 
tical uncertainties. 



previous analysis [ljj. This is explained by a 
different particle selection. In the analysis pre- 
sented here, the contribution from non-primary 
pions is larger than in the previous analysis, 
leading to smaller A and R D when using that 
procedure. However, the parameters obtained 
when applying the Bowler-Sinyukov procedure 
are almost not affected by the contribution from 
non-primary pions. 



3.4-4- Comparison of methods 

Figure 0] shows the projections of the 3- 
dimensional correlation function according to 
the Pratt-Bertsch parametrization described in 
section 21 for an azimuthally integrated anal- 
ysis. The closed symbols represent the cor- 
relation function and the open symbols the 
Coulomb corrected correlation function accord- 
ing to the standard procedure. The lines are fits 
to the data, the dashed line is the standard fit 
to the Coulomb corrected correlation function, 
and the continuous line the Bowler-Sinyukov fit 
to the uncorrected correlation function. The 
extracted parameters from both fits are the pa- 
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FIG. 6: ID correlation function for tv + tt~ com- 
pared to Standard, Bowler-Sinyukov functions and 
a theoretical calculation that includes Coulomb and 
strong interactions. 



rameters for the lowest (kx) in Fig. 

As a consistency check for the Bowler- 
Sinyukov procedure, we calculated the n + ir~ 
correlation function, dominated by Coulomb 
interaction, and compare to different calcula- 
tions. In Fig. [fj] lines indicate the standard 
(K cou i(qi nv )) and Bowler-Sinyukov ((1 — A) + 
A-Kcoui(9inv)) Coulomb functions where A was 
extracted from the fit to the 3D like-sign cor- 
relation function. This latter A is the same 
A as dilution for unlike sign pions and takes 
into account the percentage of primary pions 
through A. Clearly, the Bowler-Sinyukov func- 
tion (thick line) better reproduces the data 
(closed symbols) than the standard function 
(thin line). The small discrepancy between the 
Bowler-Sinyukov function and the data disap- 
pears when strong interaction (negligible for 
like-sign pions) is added to the Bowler-Sinyukov 
function as shown by the theoretical calculation 
[4fi| (open symbols). Between identical pions, 
there is a repulsive s-wave interaction for the 
isospin 1=2 system ^lj- However, the range 
of this interaction is estimated to be ~ 0.2 fm, 
while the characteristic separation between pi- 
ons in heavy ions collisions is ~ 5 fm. Also, 
there are no doubly charged mesonic resonances 
that could decay into same charged pions that 
would strongly interact. For these reasons, the 
strong interaction will be ignored for like sign 
particles. 



3-4-5. Coulomb interaction with the source 

The Coulomb interaction between the out- 
going charged pions and the residual positive 
charge in the source is negligible [42l |43| . This 
is confirmed by the good agreement observed 
between the parameters extracted from tt + tt + 
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FIG. 7: Momentum resolution for pions at midra- 
pidity expressed by the widths Spr/pr, Sip and 86 
as a function of p. 



and 7r 7r correlation functions as shown later 
in this paper, Fig. I14I 



3.5. Momentum resolution correction 

The limited single-particle momentum res- 
olution induces broadening of the correlation 
function and thus systematic underestimation 
of the HBT parameters. To determine the 
magnitude of this effect we need to know the 
momentum resolution for the particles under 
consideration. We estimate our single-particle 
momentum resolution by embedding simulated 
particles into real events at the TPC pixel level 
and comparing the extracted and input mo- 
menta. Figure [3 shows the RMS spreads as a 
function of \p\ in pt and angles cj> and 9, where 6 
is the angle between the momentum of the par- 
ticle and the beam axis and is the azimuthal 
angle of the particle. We see that the resolution 
inpT, given by Spx/pr (top panel of Fig-UJ) , has 
a width of about 1% for the momentum range 
under consideration. 

To account for this limited momentum res- 
olution, a correction, -Kmomcntum(<?), is applied 



to each measured correlation function: 

G(q) = — jzi rAmomcntum(gj- l 11 ) 

^ ^Plmcas ) P2meas / 

The correction factor is calculated from the 
single-particle momentum resolution as follows: 



K, 



momentum 



C((fidoal) _ B(p : 



>t(Plidoal,P2idcal) 
lidcal ,P2ideal) 



B(pismcar,P2 B mcar) 



where the ideal and smear correlation func- 
tion are formed as follows. Numerator and de- 
nominator of the ideal correlation function are 
formed by pairs of pions from different events. 
Each pair in the numerator is weighted, accord- 
ing with the Bowler- Sinyukov function, by: 



weight = (1 



(1 



A) + XK ( 



< ;oul(QWv) ^ 



where -K" C oui(<7inv) is the same factor as de- 
scribed in section 18.41 and R^ s — in the az- 
imuthally integrated analysis. If the measured 
momentum were the "real" momentum, this 
ideal correlation function would be the "real" 
correlation function. However this is not the 
case, so we calculate a smeared correlation func- 
tion for which numerator and denominator are 
also formed by pairs of pions from different 
events but their momenta have been smeared 
according to the extracted momentum resolu- 
tion. Pairs in the numerator are also weighted 
by the weight given by (|12f) . This smeared cor- 
relation function is to the ideal correlation func- 
tion, as our "measured" correlation function is 
to the "real" correlation function, which allows 
us to calculate the correction factor. 

For the weight, certain values for the HBT 
parameters (A, R , R Sl Ri and R os ) need to be 
assumed. Therefore, this procedure is iterative 
with the following steps: 

1. Fit the correlation function without mo- 
mentum resolution correction, and use 
the extracted HBT parameters for the 
first weight. 

2. Construct the momentum resolution cor- 
rected correlation function. 

3. Fit it according to Eq. (fTOft . 

4. If the extracted parameters agree with the 
parameters used to calculate the weight, 
those are the final parameters. If they 
differ from the parameters used, then use 
these latter extracted parameters for the 
new weight and go back to step 2. 
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Also, to be fully consistent, the Coulomb fac- 
tor .Kcoui(<Zinv) (where q mv is calculated from 
pairs of pions from different events) used in 
the fit to extract the HBT parameters must be 
modified to account for momentum resolution 
as follows: 



-^coul(9inv) 

= K cou i(q 



-^coul (^inv, ideal) 



inv.meas J 



: 



mv. smear 



-^coul(9mv,mcas) 



^coul(<? ] 



in v, smear 



(13) 



For this analysis, after two iterations the ex- 
tracted parameters were consistent with the in- 
put parameters. We also checked that when 
convergence is reached, the "uncorrected" HBT 
parameters matched the smeared parameters. 
The correction increases the HBT radius pa- 
rameters between 1.0% for the lowest k T bin 
[150,250] MeV/c and 2.5% for the highest bin 
[450,600] MeV/c. 



3.6. ^-dependent HBT analysis methods 

The study of HBT radii relative to the reac- 
tion plane angle was performed |44| by extend- 
ing the analysis techniques as presented in this 
section, to account for reaction plane resolu- 
tion and small instabilities in the fits. We dis- 
cuss these here. For the azimuthally-sensitive 
analysis, each of the four radii extracted from 
the Bowler-Sinyukov fit contains an implicit de- 
pendence on the azimuthal angle $ between the 
pion pair and the reaction plane. Azimuthally- 
sensitive studies of the HBT radii (i? At ($)) Q 
also must correct for finite resolution when es- 
timating the true reaction plane • Finite 
reaction plane resolution acts to decrease the 
measured amplitude of the radii oscillations, 
similar to its effect on azimuthal particle dis- 
tributions relative to the reconstructed event 
plane * 2 (i.e., elliptic flow H^)- The tech- 
nique for the resolution correction, which also 
corrects for finite $-bin width, was developed 
extensively in Ref. Q. Here we discuss briefly 
how this correction is implemented and the re- 
sulting effect on the HBT radii. 

The basic principle behind the correction 
procedure is that, for a given g-bin in the nu- 
merator A(q) and denominator B{q) of each 
correlation function, the measured contents for 
that g-bin at different <E> are modified due to 
the \& r p resolution. The true angular depen- 
dence of $ (for each <f-bin) can be extracted 
from the measured $>j by performing a Fourier 



decomposition of A(q) and B(q), which leads to 
the correction factors '9'] 



Atn(<t) = A a , n {q) 



sin(nA/2) 
nA/2 ' 



(14) 



= <„(£) (cos (n(* 2 -^ 



where a refers to both cosine and sine series, A 
is the width of each $ bin, and n is the Fourier 
component. The factors (cos (n(\E , 2 ~ x f'rp))) 
are the well-known correction factors for event 
plane resolution, obtained by extracting the 
anisotropic flow coefficients v n from the single 
particle spectrum 0,^3- The same procedure 
is used to correct the denominator B(q). 

In the present analysis, only the 2nd-order 
event plane (^2) is measured. Using Eq. 1|14|) . 
the numerator A{q) and denominator B(q) for 
each (f-bin at each measured angle <frj can be 
corrected for both the effects of angular binning 
and finite event plane resolution: 



A(q, $j 



2 £c2(A)te(glco S (2^) + 



71=1 



A e J(q)sm(2^) 



(15) 



with the correction parameter C2(A) given by 



C 2 (A) 



A 



sin(A)(cos(2(* 2 -* rp ))) p 



1. (16) 



The procedure is model independent; the quan- 
tities on the right hand side of Eq. (|15fl are all 
measured experimentally. 

For each set of $j histograms, the correc- 
tion procedure modifies both the numerators 
and denominators, and therefore the correla- 
tion functions as well. 

Figure 00 shows the squared HBT radii, ob- 
tained using Eq. (0}, as a function of $ for two 
combinations of centrality and kx- In each case, 
the oscillation amplitudes for the three trans- 
verse radii increases after the resolution correc- 
tion has been applied, while the mean shows 
little change. 

An additional technique employed in the 
azimuthally-sensitive HBT analysis is to use 
a common A parameter for each ccntrality/fcy 
bin. This step was undertaken to improve the 
quality of the fits by restricting the A parame- 
ter, under the assumption that A should have 
no implicit $ dependence. For an analysis with 
four <f> bins, this effectively reduces the number 
of free parameters per centrality/fcr bin (after 
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FIG. 8: Squared HBT radii relative to the reaction 
plane angle, without and with the reaction plane 
resolution applied, for two different centrality-fcr 
ranges. The solid lines show allowed fits to the 
individual oscillations. 



FIG. 9: Squared HBT radii relative to the reaction 
plane angle, for the case where the A parameter is 
averaged and fixed in the fit, and the case where A 
is a free fit parameter for each $. The solid lines 
show allowed fits to the individual oscillations. 



normalization) from 20 [(4 radii + A) x7V$] to 
17 [(4 radii xJV$) + A]. Since fitting all four 
correlation functions with a 17-parameter func- 
tion is arduous, we determined the average A 
parameter from the four fits and then re-fit each 
of the four correlation functions with A fixed to 
its average. 

Figure compares the fit parameters ob- 
tained with and without averaging/fixing A, for 
two centrality//cT ranges. While the individ- 
ual radii show some deviations, the resulting 
Fourier coefficients (which are represented by 
the symmetry-constrained drawn in Fig. |!5J) are 
consistent within errors for the two methods. 



3.7. Systematic uncertainties associated 
with pair cuts 

The maximum fraction of merged hits cut 
described in section 12.41 introduces a system- 
atic variation on the HBT fit parameters A, R 2 , 
R 2 S , and Rf, since it discriminates against low-g 
pairs which carry the correlation signal. This 
is a consequence of the non-Gaussianess of the 
correlation function. If it were a perfect Gaus- 
sian, this cut would not change the extracted 
parameters from the Gaussian fit, it would only 
reduce the statistics in certain bins and there- 
fore the only effect would be an increase in the 
statistical errors. 

In order to estimate this reduction we de- 
fine a range in the number of merged hits in 
which the lower limit is (i.e., no merging) and 
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FIG. 10: For ir tv (ir + Tv + ), merging occurs more 
often between tracks with \q qs\ = q q s (|g g s | 7^ 
q a q a ) than with \q q s \ =fi q q s (\q q s \ = q q s ). Note 
that exchanging the designations "1" and "2" does 
not change the sign of q q s - 



the higher limit is the value for which we con- 
sider there is too much merging. This value is 
determined from the th order Fourier coeffi- 
cients, Rq S o which is expected to be 0, Eq. JSJ. 
However, track merging introduces a deviation 
of R% a away from caused by the preferential 
merging of track pairs with correlated trans- 
verse momenta, q Q and q 8 as shown in Fig. [23 
If we calculate the components of q in the plane 
transverse to the beam as pV,i — f>T,2 where 
index 1 denotes the stiffer track and define p 
as the direction of the radius of curvature of 
the stiffer track, then if q- p is positive there is 
more merging on average. In the case of 7r~7r - 
pairs, there is a higher degree of track merg- 
ing when |(/o<7s| = Qo 

q s than when \q q s \ ^ q q s 
(top pairs). For it + tt + pairs the conditions are 
opposite (bottom pairs). 

When Rq S for tt + or ir~ analysis clearly devi- 
ates from 0, we consider that there is too much 
merging and use that value of the maximum 
fraction of merged hits as the upper limit of the 
range. We calculate the change of each HBT ra- 
dius in this range and consider that to be the 
artificial reduction due to the cut for that spe- 
cific parameter. This reduction is included as a 
systematic error in the final value. This is done 
for each centrality and each kx bin. 

As an example, Fig. ^] shows the th or- 
der (left) and 2 nd divided by th order (right) 
Fourier coefficients as a function of the maxi- 
mum fraction of merged hits allowed for the 5% 
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FIG. 11: Fourier coefficients as a function of the 
maximum fraction of merged hits for the 5% most 
central events and kr between 150 and 250 MeV/c. 



most central events and 150 < kr < 250 MeV. 
From Rq S , located in the bottom left panel, we 
determined the upper limit of merged fraction 
to be 0.2 and the corresponding variations in 
the HBT radii to be 7% for R Q , 5% for R s and 
10% for Ri. The systematic errors calculated 
according to this method are less or equal than 
10% for all radii, in all centralities and fcy bins. 



4. PION HBT AT VsNN = 200 GEV 

4.1. How Gaussian is the measured 
correlation function? 

Interferometric length scales are usually ex- 
tracted from measured correlation functions by 
fitting to a Gaussian functional form, as dis- 
cussed in [21 However, there is no reason to 
expect the measured correlation function to be 
completely Gaussian, and it is well-known that 
it seldom is. Seemingly-natural questions such 
as "how non-Gaussian is the correlation func- 
tion?" "how does the non-Gaussianness affect 
extracted length scales?" or "what is the shape 
of the correlation function?" have, unfortu- 
nately, no unique, assumption-free answers. 

Often, non-Gaussian features are simply ig- 
nored in experimental analyses. Alternatively, 
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the effect of non-Gaussianness is estimated (e.g. 
by varying the range of q- values used in the fit) 
and quoted as a systematic error on the HBT 
radii. Occasionally, some alternative functional 
form (e.g. a sum of two Gaussians, or an ex- 
ponential plus a Gaussian) is chosen ad-hoc by 
the experimenter based on the general "appear- 
ance" of the data. 

One disadvantage of this last approach is the 
difficulty in systematically comparing HBT re- 
sults obtained with different functional forms. 
Since one might hope for evidence of "new" 
physics at RHIC, it is important to place RHIC 
HBT results into the context of previously- 
established systematics. Thus this paper 
mainly focusses on Gaussian HBT radius sys- 
tematics. However, to address non- Gaussian 
issues, here we move beyond ad-hoc meth- 
ods and adopt as a standard the Edgeworth 
expansion proposed by Csorgo and collabora- 
tors |4Sl I49L l50| . Using the Bowler- Sinyukov 
Coulomb treatment, the correlation functions 
are fitted to 



FIG. 12: Projections of the 3 dimensional correla- 
tion functions and fits to Eq. (1 101 (left) and with 
the Edgeworth expansion to Eq. 11711 to 6 th order 
(right). 
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FIG. 13: HBT parameters for 0-5% most central 
events for fits to Eq. JTOJ and to Eq. JT7J to 4 th 
and to 6 th orders. Error bars reflect only statistical 
uncertainties. 



C(q Q , q s , qi) = (1 - A) + XK coui (q mv ) 
+\K coul (q inv )-e-«o R °-<>* R s-J R l X 



1+ E 



71— 4,n even 



E 

7i— 4,n even 
oc 

1+ E 



n— 4,n even 



n!(V2)" 



!(V2)» 



H n (q s R s 
H n {qiRi) , (17) 



where K i n (i = o, s, I) are fit parameters and 
H n (qiRi) are the Hermite polynomials of order 



H n (x) = (-lTe* 2 -^e-* 2 . (18) 

Only Hermite polynomials of even order are in- 
cluded in the expansion because the correlation 
function for identical particles must be invari- 
ant under (q , q s , q{) -> (-q ,-q s ,-qi)- 

At midrapidity and integrated over az- 
imuthal angle, the quantum interference term 
is factorizable into the q , g s , and qi variables. 
Therefore, it may be uniquely decomposed in 
terms of any complete set of basis functions of 
these variables. Given a sufficient number of 
terms, any basis set will do. Thus, the poten- 
tial advantage of the Edgeworth decomposition 
(Eq. I17|) is not that it is any more "model- 
independent" [4!| than, say, a Tschebyscheff 
decomposition, but that a functional expansion 
about a Gaussian shape might most econom- 
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ically describe the data. Since they are ap- 
proximately Gaussian, one hopes to capture the 
shape of measured correlation functions with 
only a few low-order terms. 

We fit our correlation functions to the form 
given by Eq. Ijl7|l for two different cases, up to 
n = 4 and up to n — 6 of the Hermite poly- 
nomials, and compare with fits to Eq. I|10|) 
(without expansion). In Fig. 1121 we show the 
fits to projections of the correlation function for 
the 0-5% most central events and kr between 
150 and 250 MeV/c, with no expansion in the 
left column and with expansion up to 6 th order 
in the right column. We observe a small im- 
provement in the fit when we include the expan- 
sion. In Fig. the extracted HBT parameters 
as a function of kr for the 0-5% most central 
events for the fits without expansion, with ex- 
pansion up to 4 th order and with expansion up 
to 6 th order are shown. In Tabled are the cor- 
responding values for the n parameters. When 
comparing the extracted parameters including 
the expansion to 6 th order to those extracted 
without the expansion, we observe that R de- 
creases by ~ 2% for all kx bins, R s changes 
between ~ -7% for the lowest kr bin [150,250] 
MeV/c and - +3% for the highest bin [450,600] 
MeV/c, and Ri decreases between ~ 18% and 
~ 8% for the lowest and highest kr bins re- 
spectively. In Table |H] are the corresponding 
X 2 /dof for those same fits. x 2 /d°f slightly im- 
proves when including the expansion up to A th 
order and does not change with the expansion 
to 6 th order. Similar trends are observed at all 
centralitics. 

We do not consider the change in HBT radii 
when including an Edgeworth expansion to rep- 
resent a systematic uncertainty when compar- 
ing to Gaussian radii traditionally discussed in 
the literature. Rather, the differences reflects 
a deviation from the Gaussian shape tradition- 
ally assumed. Furthermore, the expansion pro- 
vides a more detailed, yet still compact, charac- 
terization of the measured correlation function. 
Further theoretical development of the formal- 
ism, outside the scope of this paper, is required 
to determine whether the expansion parameters 
convey important physical information beyond 
that carried by Gaussian radius parameters. 



4.2. rriT dependence of the HBT 
parameters for most central collisions 

The HBT radius parameters measure the 
sizes of the homogeneity regions (regions emit- 
ting particles of a given momentum) |6J. Hence, 
for an expanding source, depending on the mo- 




0.5 0.6 

m T (GeV/c) 

FIG. 14: HBT parameters for 0-5% most central 
events for 7r + 7r + and -k~tx~ correlation functions. 
Error bars include statistical and systematic un- 
certainties. 



menta of the pairs of particles entering the cor- 
relation function, different parts of the source 
are measured. The size of these regions are con- 
trolled by the velocity gradients and temper- 
ature [33, loll l52j . Therefore the dependence 
of the transverse radii on transverse mass rriT 
contains dynamical information of the particle 
emitting source [H, @ . 

Figure ITU shows the HBT parameters A, R Q , 
R s , Ri and the ratio R /R s for the 0-5% most 
central events as a function of wit for 7r + 7r + 
and 7T — tt~ correlation functions. We observe 
excellent agreement between the parameters 
extracted from the positively and negatively 
charged pion analyses. The A parameter in- 
creases with rriT- This is consistent with stud- 
ies at lower energies 0, E! 13 EH > m which 
the increase was attributed to decreased con- 
tributions of pions from long-lived resonances 
at higher pr- The three HBT radii rapidly de- 
crease as a function of rtvr\ the decrease of the 
transverse radii (R and R s ) with ttit is usu- 
ally attributed to the radial flow jU IH1 E3; 
the strong decrease in Ri m ight be pro duced by 
the longitudinal flow d III El El (13 . R a falls 



16 



k T (MeV/c) 
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K;,4 
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4.92 ± 0.03 
1.25 ± 0.05 
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1.04 ± 0.06 


3.82 zb 0.04 
0.78 zb 0.06 


i?i (6 th ord.) 

Ki,4 
Ki,6 


5.01 ± 0.05 
1.32 ± 0.07 
-1.76 ± 0.29 


5.01 ± 0.04 
0.70 ± 0.07 
-2.82 ± 0.29 


4.43 ± 0.04 
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TABLE I: HBT parameters and k parameters for fits of the correlations functions without and up to 4 th 
and 6 th order of the Edgeworth expansion for the 5% most central events. The extracted radii are also 
shown in Fig. ED 



k T (MeV/c) 


No exp. 


To 4th order 


To 6th order 


150-250 


1.23 


1.09 


1.09 


250-350 


1.22 


1.05 


1.04 


350-450 


1.20 


1.02 


1.02 


450-600 


1.17 


1.01 


1.01 



TABLE II: xVdof for fits of the correlations func- 
tions without and up to 4 th and 6 th order of 
the Edgeworth expansion for the 5% most central 
events. 



steeper than R s with tot which is consistent 
with R a being more affected by radial flow [5^] . 
In contrast to many model predictions 12, 59], 
R /R s ~ 1 which indicates short emission du- 
ration in a blast wave parametrization |58l ] as 
will be discussed in next section. 

Figure 1151 compares our extracted HBT ra- 
dius parameters from 7r + 7r + and 7r~7r~ correla- 
tion functions for the 0-30% most central events 
with those obtained by the PHENIX collabora- 
tion |60| at the same beam energy and central- 
ity The same fitting procedure has been used 
in both analysis. In general, very good agree- 
ment is observed in the three radii, although 
small discrepancies are seen in R Q at small &t- 
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FIG. 15: HBT parameters from STAR and 
PHENIX at the same beam energy for the 0-30% 
most central events. Error bars include statistical 
and systematic uncertainties. 



FigurclHTlshows the HBT parameters vs. col- 
lision energy for midrapidity, lowpx 7r _ 7r _ from 
central Au+Au, Pb+Pb or Pb+Au collisions. 
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FIG. 16: Energy dependence of n~ HBT parame- 
ters for central Au+Au, Pb+Pb, and Pb+Au col- 
lisions at midrapidity and (far) ~ 0.2 GeV/c 

E3, la^, E3, E3, HE |e3, E3, Ea EE o P en 

symbols indicate that fitting was done according to 
the Bowler- Sinyukov procedure (or a similar one in 
the case of the results from Phobos). Error bars 
on NA44, NA49, CERES, PHENIX, Phobos and 
STAR results include systematic uncertainties; er- 
ror bars on other results are only statistical. 



In order to compare with our previous results 
at \Jsnn = 130 GeV, we applied similar cuts 
in our analysis as those described in 01 an d fit 
our correlation function according to the stan- 
dard procedure described in section 13.41 to ex- 
tract the HBT parameters at ^snn = 200 
GeV, closed circles at that energy in Fig. ^] 
We observe an increase of ~ 10% in the trans- 
verse radii R a and R s . In the case of R s , this 
increase could be attributed to a larger freeze- 
out volume for a larger pion multiplicity. Ri is 
consistent with our result at lower energy. The 
predicted increase by hydrodynamic models in 
the ratio R /R s as a probe of the formation 
of QGP is not observed at ^/snn = 200 GeV. 
More discussion on the lack of energy depen- 
dence of the HBT radii and its possible relation 
with the constant mean free path can be found 
in reference |65j . 
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FIG. 17: HBT parameters vs. mr for 6 different 
centralities. Error bars include statistical and sys- 
tematic uncertainties. 



We have also included in Fig. ^| the values 
for the HBT parameters at ^Jsnn = 200 GeV 
extracted when applying the cuts discussed in 
section|21and fitting the correlation function ac- 
cording to the Bowler- Sinyukov procedure (sec- 
tion 13. 4[) . open circles in the figure. This proce- 
dure is also used by the CERES collaboration. 
The smaller A, R Q , and R s can be explained by 
the different cuts as already discussed in sec- 
tion The larger value for R a /R s is due to 
the improved procedure of taking Coulomb in- 
teraction into account in the Bowler- Sinyukov 
procedure, section EOl 



4.3. Centrality dependence of the mr 
dependence 

We observe excellent agreement between the 
results for positively and negatively charged 
pion correlation functions for the most central 
collisions shown in the previous section. There- 
fore, we add the numerators and denominators 
of the correlation functions for positive and neg- 
ative pions in order to improve statistics; all 
the results shown in the rest of this section cor- 
respond to these added correlation functions. 
The centrality dependence of the source param- 
eters is presented in Fig.ll7lwhere the HBT pa- 
rameters are shown as a function of mx for 6 
different centralities. The A parameter slightly 
increases with decreasing centrality. The three 
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FIG. 18: Fourier coefficients of azimuthal oscilla- 
tions of HBT radii vs. number of participating nu- 
cleons, for tt + and tt~ pairs separately (0.25 < kr < 
0.35 GeV/c). Left panels: means (0 th -order FC) of 
oscillations; right panels: relative amplitudes (see 
text for details). Larger participant numbers cor- 
respond to more central collisions. 



radii increase with increasing centrality and Ri 
varies similar to R a and R s . For R Q and R s this 
increase may be attributed to the initial geo- 
metrical overlap of the two nuclei. R /R s ~ 1, 
for all centralities. 



4.4. Azimuthally-sensitive HBT 

The results presented in Q were for 7r+ and 
7r~ correlation functions combined before fit- 
ting. Figure ED shows the consistency of the 
Fourier coefficients obtained with Eq. (JSJ as a 
function of number of participants. 

In section [3] we noted that the th order 
Fourier coefficients correspond to the extracted 
HBT radii in an azimuthally integrated anal- 
ysis. This is confirmed in Fig. EH that shows 
the excellent agreement between them. The 
azimuthally integrated [traditional) HBT radii 
(closed symbols) agree within 1/10 fm with the 
th -order Fourier coefficients (open symbols) 
from the azimuthally-sensitive HBT analysis. 



5. DISCUSSION OF RESULTS 

As already mentioned in section FOl the de- 
pendence of the HBT radii on the transverse 
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FIG. 19: Comparison between the HBT radii 
obtained from and azimuthally integrated (tradi- 
tional) HBT analysis and the th -order Fourier co- 
efficients from an azimuthally-sensitive HBT anal- 
ysis. Error bars include only statistical uncertain- 
ties. 
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FIG. 20: HBT radius parameters for 6 differ- 
ent centralities. The lines indicate power-law fits 
(Ri(rriT) = R'i ■ (mT/m 7T )~ c " i ) to each parameter 
for each centrality. 



mass contains dynamical information about the 
particle emitting source. In order to extract 
this information, we fit the mr dependence of 
the HBT radii for each centrality from Fig. IT7I 
using a simple power-law fit: ^(tot) = R'i • 
(mr /m 7T )~ ai (solid lines in Fig. I20fl . Figure 
[?T1 shows the extracted fit parameters for the 
three HBT radii, R' in the top panel and a in 
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radius parameters (lines in Fig. I2UI . 



the lower panel, as a function of the number of 
participants, where A^ par ticipants has been calcu- 
lated from a Glauber model described in [6(|. 
-^participants increases with the centrality of the 
collision. R' decreases with decreasing number 
of participants, which is consistent with the de- 
creasing initial source size, a is approximately 
constant for Ri which would indicate that the 
longitudinal flow is similar for all centralitics. 
However, for the transverse radii R and R s , a 
seems to decrease for the most peripheral col- 
lisions which could be an indication of a small 
reduction of transverse flow and/or an increase 
of temperature for those most peripheral colli- 
sions. This is consistent with the values for flow 
and temperature extracted from blast wave fits 
to pion, kaon, and proton transverse momen- 
tum spectra [63, as well as to HBT as will be 
discussed in next section. The drop of a with 
decreasing number of participants is faster in 
R than in R s which could again indicate that 
R might be more affected by radial flow [Hsj . 



5.1. Blast wave parametrization 

Hydrodynamic calculations that successfully 
reproduce transverse momentum spectra and 
elliptic flow, fail to reproduce the HBT param- 
eters [5^. In most cases, these calculations un- 
derestimate R s and overestimate R and i?;. 



Since R s only probes the spatial extent of the 
source while R and Ri are also sensitive to the 
system lifetime and the duration of the parti- 
cle emission they may be underestimating 
the system size and overestimating its evolution 
time and emission duration. We fit our data 
with a blast wave parametrization designed to 
describe the kinetic freeze-out configuration. In 
this section we will discuss the extracted pa- 
rameters and their physical implications. 

This blast wave parametrization jE^ assumes 
that the system is contained within an infinitely 
long cylinder along the beam line and requires 
longitudinal boost invariant flow. It will be 
shown that this latter assumption is not nec- 
essarily correct. It also assumes uniform parti- 
cle density. The single set of free parameters in 
this parametrization is: the kinetic freeze-out 
temperature (T); the maximum flow rapidity 
(p = f(po + p a cos(2r/>)), for an azimuthally in- 
tegrated analysis p a = 0); the radii (R for the 
azimuthally integrated analysis and R x , R y for 
the azimuthally sensitive analysis) of the cylin- 
drical system; the system longitudinal proper 
time (t = \/t 2 — z 2 )\ and the emission dura- 
tion (At). 

We use this parametrization to fit the az- 
imuthally integrated pion HBT radii, as well 
as the azimuthally sensitive pion HBT radii. In 
both fits, T and po are fixed to those extracted 
from a blast wave fit to pion, kaon and proton 
transverse momentum spectra |67j and 112 |68|. 
By doing this, the azimuthally integrated and 
azimuthally sensitive radii are fitted with the 
same temperature and po, and the edge source 
radii can be compared directly. Also, a 5% er- 
ror was added to all HBT radii before the fit in 
order to reflect the blast wave systematic errors 
described in J§^. In the fit, the transverse flow 
rapidity linearly increases from zero at the cen- 
ter to a maximum value at the edge of the sys- 
tem. The best fit parameters are summarized in 
Table ITTT1 for the azimuthally integrated analy- 
sis, and Table llVl for the azimuthally sensitive 
analysis. 

Most of the parameters, as well as their evo- 
lution with centrality, agree with similar stud- 
ies. Temperature decreases with increasing cen- 
trality and the average transverse flow velocity 
((Pt) = J &rct&nh(p ^)rdr / J rdr) increases 
with increasing centrality. Both results are con- 
sistent with those extracted from fits to spec- 
tra only |67j and reflect increased rescattering, 
expansion and system evolution time with in- 
creasing centrality. 

Figure [23 shows the R parameter extracted 
from the blast wave fit as they are in Table ITTT1 
Also shown in that plot is i? ge0 m calculated as- 
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Centrality (%) 


T (MeV) 


Po 


R (fm) 


r (fm/c) 


At (fm/c) 


X 2 /dof 


0-5 


97 ± 2 


1.03 ± 0.01 


13.3 ± 0.2 


9.0 ± 0.3 


2.83 ± 0.19 


3.13/9 


5-10 


98 ± 2 


1.00 ± 0.01 


12.6 ± 0.2 


8.7 ± 0.2 


2.45 ± 0.17 


2.71/9 


10-20 


98 ± 3 


0.98 ± 0.01 


11.5 ± 0.2 


8.1 ± 0.2 


2.35 ± 0.16 


2.61/9 


20-30 


100 ± 2 


0.94 ± 0.01 


10.5 ± 0.1 


7.2 ± 0.1 


2.10 ± 0.09 


0.99/9 


30-50 


108 ± 2 


0.86 ± 0.01 


8.8 ± 0.1 


5.9 ± 0.1 


1.74 ± 0.12 


2.13/9 


50-80 


113 ± 2 


0.74 ± 0.01 


6.5 ± 0.1 


4.0 ± 0.2 


1.73 ± 0.10 


1.12/9 



TABLE III: Extracted parameters from a blast wave fit to azimuthally integrated pion HBT radii, with T 
and po fixed from fits to pion, kaon, and proton transverse momentum spectra and vi. 



Cent 


(%) 


T (MeV) 


Po 


P2 


Ra, (fm) 


R y (fm) 


t (fm/c) 


At (fm/c) 


X 2 /dof 


0- 


-5 


97 ± 2 


1.03 ± 0.01 


0.03 ± 0.002 


12.9 ± 0.1 


13.4 ± 0.1 


8.9 ± 0.2 


3.16 ± 0.11 


106.8/63 


5- 


10 


98 ± 2 


1.00 ± 0.01 


0.04 ± 0.002 


12.1 ± 0.1 


12.9 ± 0.1 


8.2 ± 0.2 


2.73 ± 0.12 


103.2/63 


10- 


-20 


98 ± 3 


0.98 ± 0.01 


0.05 ± 0.002 


10.9 ± 0.1 


11.9 ± 0.1 


7.8 ± 0.2 


2.59 ± 0.10 


131.06/63 


20- 


-30 


100 ± 1 


0.94 ± 0.01 


0.07 ± 0.002 


9.7 ± 0.1 


11.0 ± 0.1 


6.9 ± 0.1 


2.29 ± 0.10 


87.2/63 


30- 


-80 


112 ± 2 


0.82 ± 0.05 


0.10 ± 0.005 


7.4 ± 0.1 


8.7 ± 0.1 


5.1 ± 0.2 


1.94 ± 0.14 


189.4/63 



TABLE IV: Extracted parameters from a blast wave fit to azimuthally sensitive pion HBT radii, with T 
and po fixed from fits to pion, kaon, and proton transverse momentum spectra and vi. 
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FIG. 22: Extracted freeze-out source radius ex- 
tracted from a blast wave fit; source radius i? gC0 m 
from fits to R a (lines in Fig. I23II : and 2 • R s for 
the lowest kr bin as a function of number of par- 
ticipants. R from blast wave contains only uncer- 
tainties from fit; R ga0 m error bars contain system- 
atic uncertainties from the input parameters; 2 • R s 
contains the statistical and systematic uncertain- 
ties from R s . 



suming a transverse expanding, longitudinally 
boost-invariant source, and a Gaussian trans- 
verse density profile, by fitting the tot depen- 
dence of R s to 52]: 



R s {m T ) 




R 2 

''geom 



(19) 



where T is the freeze-out temperature and po is 
the surface transverse rapidity. Figurel2"51shows 
such fits to R s for each centrality with T and 
Po extracted from blast wave fits to pion, kaon, 
and proton transverse momentum spectra (T 




0.5 0.6 

m T (GeV/c) 



FIG. 23: HBT parameter R 3 . Lines represent the 

fits R s (m T ) = ] jRlom/[l+p 2 o^ + ^)]- 



= 90 MeV, po — 1.20 for the most central colli- 
sions and T = 120 MeV, p = 0.82 for the most 
peripheral bins) [63]. Fig. [23 shows good agree- 
ment between these two extracted radii that in- 
crease from ~ 5 fm for the most peripheral col- 
lisions to ~ 13 fm for the most central collisions 
following the growth of the system initial size. 
The differences may be explained by the poor 
quality of the fits to R s as seen in Fig. 1231 

As already mentioned, R s carries only spa- 
tial information about the source |32l l33j] . In 
the special case of vanishing space-momentum 
correlations (no transverse flow or T — ► oo), the 
source spatial distribution may be modelled by 
a uniformly filled disk of radius R, which in 
this case is exactly two times R s , the RMS of 
the distribution along a specific direction. In 
Fig. |2Uwe have included 2 • R s for our lowest kx 
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Centrality 


Rx, initial (fill) 


Ry, initial (fill) 


0-5% 


5.70 ± 0.01 


5.86 ± 0.01 


5-10% 


5.28 ± 0.01 


5.72 ± 0.01 


10-20% 


4.74 ± 0.01 


5.50 ± 0.01 


20-30% 


4.14 ± 0.01 


5.12 ± 0.01 


30-50% 


3.58 ± 0.01 


4.70 ± 0.01 


50-80% 


2.84 ± 0.01 


4.02 ± 0.01 


30-80% 


3.48 ± 0.01 


4.60 ± 0.01 



TABLE V: Initial in-plane (.fi^, initial) and out-of- 
plane (R Vi initial) radii for 7 different centrality bins. 



bin, kr = [150,250] MeV/c, in order to compare 
it with the extracted source radii. We observe 
the effect of space-momentum correlations that 
reduce the size of the regions of homogeneity in 
the results from the most central collisions for 
which 2 • R s is smaller than the extracted radii 
from the fit. 

For an azimuthally asymmetric collision, the 
initial source has an elliptic shape with the 
larger axis perpendicular to the reaction plane 
(out-of-plane) and the shorter axis in the re- 
action plane (in-plane). In order to calculate 
the radii of the initial source in the x (in- 
plane) and y (out-of-plane) direction we first 
get the initial distribution of particles in the 
almond shaped initial overlap from a Monte 
Carlo Glauber model calculation as described 
in [6(J. The in-plane (R x , initial) and out-of- 
plane (R y , initial) initial radii are calculated as 
the radii of the region that contains 95% of the 
particles. The values for the initial in-plane 
and out-of-plane edge radii are shown in Ta- 
ble[Vl The azimuthally integrated initial radius 
(-^initial) can be calculated from those two radii 
as: 



R 



.initial 



initial 



R 2 

y, initial 



(20) 



Figure |2] (bottom panel) shows i?/i?i n itiai 
vs. number of participants for in-plane, out- 
of-plane and azimuthally integrated directions. 
The final source radii are those extracted from 
the blast wave parametrization. i?/i?i n itiai is 
the relative expansion of the source which is 
stronger in-plane than out-of-plane for the most 
peripheral collisions, and it is similar in both 
directions for the most central collisions. The 
azimuthally integrated radius indicates a strong 
relative expansion of the source for central col- 
lisions. This expansion seems to be very sim- 
ilar for all centralities, decreasing just for the 
most peripheral cases. Figure (top panel) 
shows the overall expansion of the source given 
by R — i?initiai vs. number of participants. 
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FIG. 24: R-Rmitial (top panel) and i?/i? in itiai (bot- 
tom panel) for the azimuthally integrated analysis 
and in the x (in-plane) and y (out-of-plane) direc- 
tions for the azimuthally-sensitive case vs. number 
of participants. 
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FIG. 25: R — -Rinitiai for the azimuthally integrated 
analysis and in the x (in-plane) and y (out-of- 
plane) directions for the azimuthally-sensitive case 
vs. (dN/dy)/R\ 

nitial - 



While the absolute expansion (R— -Rinitiai) in- 
creases steadily going to more central collision, 
the relative expansion saturates when the num- 
ber of participant reaches 150. Furthermore, 
both absolute and relative expansions differ sig- 
nificantly in peripheral events when comparing 
the in-plane and out-of-plane directions, while 
they are similar for the most central. This is ex- 
pected arguing that the expansion, or in other 
word flow, is driven by particle reinteractions 
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following the initial pressure gradients which 
in turn follow the initial energy density gradi- 
ents. As the centrality increases the difference 
between the initial energy density gradient in- 
plane and out-of-plane diminishes, which brings 
the expansions in-plane and out-of-plane closer 
together. 

The question is then what drives the trans- 
verse expansion. The difference between the 
in-plane and out-of-plane expansion at a given 
centrality shows that the initial energy density 
gradient matters. The initial energy density 
gradients are responsible for establishing the 
initial expansion velocity but the spatial expan- 
sion will also depend on for how long the system 
expands. The system lifetime is likely to de- 
pend on the initial energy density, which may 
be gauged by dividing the particle multiplic- 
ity (dN/dy) by the initial area estimated in the 
Glauber framework as done in |6fij . Following 
this idea, we investigate how the transverse ex- 
pansion evolves while varying centrality, which 
affects both the average energy density and the 
energy density gradient. We find that the trans- 
verse expansion scales with (dN / dy) / Ri u iti a \ as 
shown in Fig. [23 This figure shows R — i?i n itiai 
vs. (dN/dy) /RinMt,], where dN/dy is for pi- 
ons as reported in |67| and -Rinitiai is the corre- 
sponding in-plane, out-of-plane or azimuthally 
integrated initial radius described above. This 
quantity scales neither as a gradient nor as an 
energy density but it appears to contain the rel- 
evant parameters that drive the transverse ex- 
pansion. We observe a clear scaling for R x and 
R y as well as for the azimuthally integrated ra- 
dius R, with (dN/dy) /Ri n iti a i- For the same 
collisions, the in-plane expansion corresponds 
to a higher value of (dN / dy) / -R; n itiai than the 
corresponding out-of-plane expansion. 

The good fit to the data obtained with the 
blast wave parametrization, consistent with ex- 
pansion, and the comparison in different ways 
of the initial and final sizes of the source clearly 
indicate that the results can be interpreted 
in terms of collective expansion that could be 
driven by the initial pressure gradient. How- 
ever, the time scales extracted from the fit seem 
to be very small, smaller than the values pre- 
dicted by hydrodynamic models. 

From the dependence of Ri on ttlt shown in 
Fig. 1261 and assuming boost-invariant longitu- 
dinal flow, we can extract information about 
the evolution time-scale of the source, or proper 
time of freeze-out, by fitting it to a formula first 
suggested by Sinyukov and collaborators 0, |5(j 
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FIG. 26: Longitudinal HBT radius Ri. Lines rep- 
resent the fits Ri = T\ 
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FIG. 27: Evolution time r vs. number of partici- 
pants as extracted from a fit to Ri, lines in Figure 
1261 (triangles) , and from a blast wave fit to HBT 
parameters and spectra (circles). 



and then improved by others |58j : 



Ri 



I T K 2 (m T /T) 
rriT Ki(rriT /T) 



(21) 



where T is the freeze-out temperature and K\ 
and K2 are the modified Bessel functions of or- 
der 1 and 2. This expression for Ri also as- 
sumes vanishing transverse flow and instanta- 
neous freeze-out in proper time (i.e., At = 0). 
The first assumption is approximatively justi- 
fied by the small dependence of Ri on po in 
full calculation 58] . The second approximation 
is justified by the small Ar from blast wave 
fits, Table II I II Figure [25] also shows the fits 
to Ri (lines) using temperatures, T, consistent 
with spectra as for the fit to R s . The extracted 
values for the evolution time r are shown in 
Fig. 1271 The evolution time increases with cen- 
trality from t m 4 fm/c for the most periph- 
eral events to r Ri 9 fm/c for the most central 
events. In the same plot, the extracted evo- 



2.3 



In plane (R„) 
Out of plane (R y ) 
Azim. Integrated (R) 




^3.5 
o 

I 3 



<2.5 
2 
1.5 
1 

0.5 




-c(fm) 



o 



50 100 150 200 250 



300 N 350 

participants 



400 



FIG. 28: 7? — -Rmitiai for the azimuthally integrated 
analysis and for the in-plane and out-of-plane di- 
rections vs. /?T,max • t. The line is a "y — x" line. 



FIG. 29: Emission duration time At vs. number of 
participants as extracted using a blast fit to HBT 
parameters and spectra. 



lution time from the blast wave fit is shown. 
Good agreement is observed between the two 
extracted proper times for all centralities. They 
are surprisingly small as compared with hydro- 
dynamical calculations that predict a freeze-out 
time of ~ 15 fm/c in central collisions. These 
hydrodynamical calculations may over-predict 
the system lifetime or the assumption on which 
the extraction of r is based in the blast wave 
parametrization, longitudinal boost invariant 
expansion, might not be completely justified. 

As a check for the consistency of the evo- 
lution time extracted from the blast wave fit, 
Fig. EHl shows the final source radius as ex- 
tracted from the blast wave fit minus the ini- 
tial source size vs. /Sr.max ■ t. This /Sr.max is 
the maximum flow velocity and is expected to 
be the velocity at the edge of the expanding 
source at kinetic freeze-out. It has been calcu- 
lated from the po and p a blast wave parameters 
as /?T,max = tanh(p + p a cos(2f/>)). <f> is in- 
plane and 7r/2 out-of-plane £3, and p a is for 
the azimuthally integrated analysis and is given 
in Table llVl for the azimuthally sensitive anal- 
ysis. The evolution time, r, is the blast wave 
parameter shown in Fig. 1271 and Table ITTT1 The 
systematic errors in /3T,max ■ t come from the 
finite size bin in centrality. If the extracted ra- 
dius and proper-time are right, the initial and 
final edge radii should be related by the rela- 
tion i?fi nal < ^initial + /3T,max ' T SO that 

the points in the figure should all be clearly be- 
low the solid line (Ar.max -t = R- initial)- 
Since most points are above the line, a possible 
explanation is that r is not properly calculated 
within the blast wave parametrization. A larger 
r would move the points below the line. 

Figure TM shows the emission duration time, 
At as a function of number of participants. At 
increases with increasing centrality up to ~ 3 
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FIG. 30: Final source eccentricity (eanai) as cal- 
culated from the Fourier coefficients (2i2 Si2 /^so) 
and from the final in-plane and out-of-plane radii 
((Ry — R x )/(Ry + Rx)) vs - initial eccentricity 
(einitiai). The most peripheral collisions corre- 
spond to the largest eccentricity. The line indicates 
£flnai = einitiai- Systematic errors of 30%, based on 
sensitivity to model parameters |58||. are assigned 
to £fl na i extracted from the Fourier coefficients. 



fm/c. It is relatively small for all centralities, 
however it has increased with respect to the val- 
ues extracted from our analysis at \J snn = 130 
GeV |^ due to the improved procedure of tak- 
ing Coulomb interaction into account and the 
consequent increase in R Q . 

The freeze-out shape of the source in non- 
central collisions and its relation to the spatial 
anisotropy of the collision's initial overlap re- 
gion give us another hint about the system life- 
time. The initial anisotropic collision geome- 
try generates greater transverse pressure gradi- 
ents in the reaction plane than perpendicular 
to it. This leads to a preferential in-plane ex- 
pansion 0, [2(], El which diminishes the ini- 
tial anisotropy. A long-r source would be less 
out-of-plane extended and perhaps in-plane ex- 
tended. The eccentricity of the initial over- 
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lap region has been calculated from the initial 
RMS of the distribution of particles, as given 
by the Monte Carlo Glauber calculation, in the 
in-plane and out-of-plane directions as: 

(nRMS \2 /dRMS \2 

V y. initial/ V initial/ fOO^ 

£ = ( nRMS \2 7 ( pRMS \2 ' V - 1 
VV initial-' v Jt a:, initial^ 

Figure |201 shows the relation between initial 
and final eccentricities, with more peripheral 
collisions showing a larger final anisotropy. The 
final source eccentricity has been calculated 
from the Fourier coefficients {2R 2 2 /R 2 .o)j as 
well as from the final in-plane and out-of-plane 
radii {{R 2 y - R 2 X ) / {R 2 y + R 2 X )) extracted from the 
blast wave fit to azimuthally sensitive HBT and 
spectra described above. The source at freeze- 
out remains out-of-plane extended indicating 
that the outward pressure and/or expansion 
time was not sufficient to quench or reverse the 
initial spatial anisotropy. The large elliptic flow 
and small HBT radii observed at RHIC energies 
might favor a large pressure build-up in a short- 
lived system compared to hydrodynamic cal- 
culations. Also, out-of-plane freeze-out shapes 
tend to disfavor a long-lived hadronic rescat- 
tering phase following hydrodynamic expansion 
|23j . This short hadronic phase is consistent 
with a short emission time. 



6. CONCLUSION 

We have presented a detailed description of 
a systematic HBT analysis in Au+Au collisions 
at yfs~N~N = 200 GeV. We have analyzed the 



Gaussianess of the correlation function and con- 
clude that there is a deviation from a pure 
Gaussian, however it is not clear how this af- 
fects the HBT parameters, extracted assuming 
a Gaussian, that will be compared to models. 
We have studied the centrality dependence of 
the kr dependence of the HBT parameters and 
extracted geometrical and dynamical informa- 
tion on the source at freeze-out. We conclude 
that there is a significant expansion in Au+Au 
collisions, and that the relative expansion does 
not significantly depend on centrality. The sys- 
tem expands by a factor of at least 2.0 for most 
centralities. This is well established by HBT. 
The initial pressure gradient seems to be driv- 
ing the expansion. The extracted time scales 
from a blast wave fit are small. The blast wave 
evolution time r is small as compared with hy- 
drodynamical calculations which could suggest 
that the longitudinal boost invariant assump- 
tion has only limited validity. 
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